using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.QuadPack
{
    /// <summary>
    /// DQC25O - Integration rules for functions with cos or sin factor.
    /// </summary>
    [Serializable]
    public static class Dqc25oClass
    {
        /// <summary>
        /// Integration rules for functions with cos or sin factor.
        /// </summary>
        /// <param name="f">The function.</param>
        /// <param name="a">The lower limit of integration.</param>
        /// <param name="b">The upper limit of integration.</param>
        /// <param name="omega">Double precision parameter in the weight function.</param>
        /// <param name="sincos">Indicates which weight function is to be used: = 1 - W(X) = COS(OMEGA*X); = 2 - W(X) = SIN(OMEGA*X).</param>
        /// <param name="nrmom">The length of interval (a,b) is equal to the length of the original integration interval divided by 2**NRMOM.</param>
        /// <param name="maxp1">Gives an upper bound on the number of chebyshev moments which can be stored, i.e. for the intervals of lengths ABS(BB-AA)*2**(-L), L = 0,1,2, ..., MAXP1-2.</param>
        /// <param name="ksave">Key which is one when the moments for the current interval have been computed.</param>
        /// <param name="abserr">Estimate of the modulus of the absolute error, which should not exceed abs(i-result).</param>
        /// <param name="neval">The needed function evalations.</param>
        /// <param name="resabs">Approximation to the integral j.</param>
        /// <param name="resasc">Approximation to the integal of abs(f-i/(b-a)) over (a,b).</param>
        /// <param name="momcom">The momcom parameter.</param>
        /// <param name="chebmo">The chebmo parameter.</param>
        /// <returns>Approximation to the integral.</returns>
        public static double Dqc25o(HardRealFunction f, double a, double b, double omega, int sincos, int nrmom,
                                    int maxp1, int ksave, ref double abserr, ref int neval, ref double resabs,
                                    ref double resasc, ref int momcom, double[,] chebmo)
        {
            double[] x = {
                             0.99144486137381041114, 0.96592582628906828675, 0.92387953251128675613,
                             0.86602540378443864676
                             , 0.79335334029123516458, 0.70710678118654752440, 0.60876142900872063942,
                             0.50000000000000000000, 0.38268343236508977173, 0.25881904510252076235,
                             0.13052619222005159155
                         };

            double ac;
            double an;
            double an2;
            double @as;
            double asap;
            double ass;
            double centr;
            double conc;
            double cons;
            double cospar;
            double estc;
            double ests;
            double hlgth;
            int mm1 = 0;
            double parint;
            double par2;
            double par22;
            double p2 = 0;
            double p3 = 0;
            double p4 = 0;
            double resc12;
            double resc24;
            double ress12;
            double ress24;
            double result;
            double sinpar;
            double[] cheb12 = new double[13];
            double[] cheb24 = new double[25];
            double[] c = new double[28];
            double[] d = new double[28];
            double[] d1 = new double[28];
            double[] d2 = new double[28];
            double[] d3 = new double[28];
            double[] fval = new double[25];
            double[] v = new double[28];
            int i;
            int isym;
            int j;
            int k;
            int m;
            int noequ;
            int noeq1;

            centr = 0.5 * (b + a);
            hlgth = 0.5 * (b - a);
            parint = omega * hlgth;

            // Compute the integral using the 15-point Gauss-Kronrod formula
            // * if the value of the parameter in the integrand is small or
            // * is less than (bb-aa)/2^(maxp1-2), where (aa,bb) is the original
            // * integration interval.
            // 
            if (Math.Abs(parint) > 2.0)
            {
                goto _10;
            }
            result = Dqk15wClass.G_K15W(f, DqwgtClass.Dqwgto, omega, p2, p3, p4, sincos, a, b, ref abserr, ref resabs,
                                        ref resasc);
            neval = 15;
            goto _190;

            // Compute the integral using the generalized Clenshaw-Curtis method. 
            _10:
            conc = hlgth * Math.Cos(centr * omega);
            cons = hlgth * Math.Sin(centr * omega);
            resasc = Double.MaxValue;
            neval = 25;

            // Check whether the Chebyshev moments for this interval have
            // * already been computed.
            // 
            if ((nrmom < momcom) || (ksave == 1))
            {
                goto _140;
            }

            // Compute a new set of Chebyshev moments. 
            m = momcom + 1;
            //** Add variable mm1 to ease transliteration from FORTRAN array
            // *** indexing to C indexing.
            // **
            mm1 = m - 1;
            par2 = parint * parint;
            par22 = par2 + 2.0;
            sinpar = Math.Sin(parint);
            cospar = Math.Cos(parint);

            // Compute the Chebyshev moments with respect to cosine. 
            v[0] = 2.0 * sinpar / parint;
            v[1] = (8.0 * cospar + (par2 + par2 - 8.0) * sinpar / parint) / par2;
            v[2] = (32.0 * (par2 - 12.0) * cospar + (2.0 * ((par2 - 80.0) * par2 + 192.0) * sinpar) / parint) /
                   (par2 * par2);
            ac = 8.0 * cospar;
            @as = 24.0 * parint * sinpar;
            if (Math.Abs(parint) > 24.0)
            {
                goto _70;
            }

            // Compute the Chebyshev moments as the solutions of a boundary
            // * value problem with 1 initial value (v[2]) and 1 end value
            // * (computed using an asymptotic formula).
            // 

            noequ = QuadPackConstants.NMAC - 3;
            noeq1 = noequ - 1;
            an = 6.0;
            for (k = 0; k <= noeq1; k++)
            {
                an2 = an * an;
                d[k] = -2.0 * (an2 - 4.0) * (par22 - an2 - an2);
                d2[k] = (an - 1.0) * (an - 2.0) * par2;
                d1[k] = (an + 3.0) * (an + 4.0) * par2;
                v[k + 3] = @as - (an2 - 4.0) * ac;
                an += 2.0;
            }
            an2 = an * an;
            d[noequ] = -2.0 * (an2 - 4.0) * (par22 - an2 - an2);
            v[noequ + 3] = @as - (an2 - 4.0) * ac;
            v[3] -= (56.0 * par2 * v[2]);
            ass = parint * sinpar;
            asap = (((((210.0 * par2 - 1.0) * cospar - (105.0 * par2 - 63.0) * ass) / an2 - (1.0 - 15.0 * par2) * cospar +
                      15.0 * ass) / an2 - cospar + 3.0 * ass) / an2 - cospar) / an2;
            v[noequ + 3] -= (2.0 * asap * par2 * (an - 1.0) * (an - 2.0));
            // Solve the tridiagonal system by means of Gaussian elimination
            // * with partial pivoting.
            // 
            for (i = 0; i <= noequ; i++)
            {
                d3[i] = 0.0;
            }
            d2[noequ] = 0.0;
            for (i = 0; i <= noeq1; i++)
            {
                if (Math.Abs(d1[i]) <= Math.Abs(d[i]))
                {
                    goto _40;
                }
                an = d1[i];
                d1[i] = d[i];
                d[i] = an;
                an = d2[i];
                d2[i] = d[i + 1];
                d[i + 1] = an;
                d3[i] = d2[i + 1];
                d2[i + 1] = 0.0;
                an = v[i + 4];
                v[i + 4] = v[i + 3];
                v[i + 3] = an;
                _40:
                d[i + 1] -= (d2[i] * d1[i] / d[i]);
                d2[i + 1] -= (d3[i] * d1[i] / d[i]);
                v[i + 4] -= (v[i + 3] * d1[i] / d[i]);
            }
            v[noequ + 3] /= d[noequ];
            v[noequ + 2] = (v[noequ + 2] - d2[noeq1] * v[noequ + 3]) / d[noeq1];
            for (i = 1; i <= noeq1; i++)
            {
                k = noequ - i - 1;
                v[k + 3] = (v[k + 3] - d3[k] * v[k + 5] - d2[k] * v[k + 4]) / d[k];
            }
            goto _90;

            // Compute the Chebyshev moments by means of forward recursion. 
            _70:
            an = 4.0;
            for (i = 3; i < 13; i++)
            {
                an2 = an * an;
                v[i] = ((an2 - 4.0) * (2.0 * (par22 - an2 - an2) * v[i - 1] - ac) + @as -
                        par2 * (an + 1.0) * (an + 2.0) * v[i - 2]) / (par2 * (an - 1.0) * (an - 2.0));
                an += 2.0;
            }
            _90:
            for (j = 0; j < 13; j++)
            {
                chebmo[mm1, 2 * j] = v[j];
            }

            // Compute the Chebyshev moments with respect to sine. 
            v[0] = 2.0 * (sinpar - parint * cospar) / par2;
            v[1] = (18.0 - 48.0 / par2) * sinpar / par2 + (-2.0 + 48.0 / par2) * cospar / parint;
            ac = -24.0 * parint * cospar;
            @as = -8.0 * sinpar;
            chebmo[mm1, 1] = v[0];
            chebmo[mm1, 3] = v[1];
            if (Math.Abs(parint) > 24.0)
            {
                goto _120;
            }
            for (k = 2; k < 12; k++)
            {
                an = k + 1;
                chebmo[mm1, 2 * k + 1] = -sinpar / (an * (2.0 * an - 2.0)) -
                                         0.25 * parint * (v[k + 1] / an - v[k] / (an - 1.0));
            }
            goto _140;

            // Compute the Chebyshev moments by means of forward recursion. 
            _120:
            an = 3.0;
            for (i = 2; i < 12; i++)
            {
                an2 = an * an;
                v[i] = ((an2 - 4.0) * (2.0 * (par22 - an2 - an2) * v[i - 1] + @as) + ac -
                        par2 * (an + 1.0) * (an + 2.0) * v[i - 2]) / (par2 * (an - 1.0) * (an - 2.0));
                an += 2.0;
                chebmo[mm1, 2 * i + 1] = v[i];
            }
            _140:
            if (nrmom < momcom)
            {
                m = nrmom + 1;
                mm1 = m - 1;
            }
            if ((momcom < (maxp1 - 1)) && (nrmom >= (momcom)))
            {
                (momcom)++;
            }

            // Compute the coefficients of the Chebyshev expansions of degrees
            // * 12 and 24 of the function f.
            // 
            fval[0] = 0.5 * f.SolveAt(centr + hlgth);
            fval[12] = f.SolveAt(centr);
            fval[24] = 0.5 * f.SolveAt(centr - hlgth);
            for (i = 1; i < 12; i++)
            {
                isym = 24 - i;
                fval[i] = f.SolveAt(hlgth * x[i - 1] + centr);
                fval[isym] = f.SolveAt(centr - hlgth * x[i - 1]);
            }

            DqchebClass.Dqcheb(x, fval, cheb12, cheb24);

            // Compute the integral and error estimates. 
            resc12 = cheb12[12] * chebmo[mm1, 12];
            ress12 = 0.0;
            estc = Math.Abs(cheb24[24] * chebmo[mm1, 24]) + Math.Abs((cheb12[12] - cheb24[12]) * chebmo[mm1, 12]);
            ests = 0.0;
            k = 10;
            for (j = 0; j < 6; j++)
            {
                resc12 += (cheb12[k] * chebmo[mm1, k]);
                ress12 += (cheb12[k + 1] * chebmo[mm1, k + 1]);
                estc += Math.Abs((cheb12[k] - cheb24[k]) * chebmo[mm1, k]);
                ests += Math.Abs((cheb12[k + 1] - cheb24[k + 1]) * chebmo[mm1, k + 1]);
                k -= 2;
            }
            resc24 = cheb24[24] * chebmo[mm1, 24];
            ress24 = 0.0;
            resabs = Math.Abs(cheb24[24]);
            k = 22;
            for (j = 0; j < 12; j++)
            {
                resc24 += (cheb24[k] * chebmo[mm1, k]);
                ress24 += (cheb24[k + 1] * chebmo[mm1, k + 1]);
                resabs += (Math.Abs(cheb24[k]) + Math.Abs(cheb24[k + 1]));
                if (j <= 4)
                {
                    estc += (Math.Abs(cheb24[k] * chebmo[mm1, k]));
                    ests += (Math.Abs(cheb24[k + 1] * chebmo[mm1, k + 1]));
                }
                k -= 2;
            }
            resabs *= Math.Abs(hlgth);
            if (sincos == QuadPackConstants.SINE)
            {
                goto _180;
            }
            result = conc * resc24 - cons * ress24;
            abserr = Math.Abs(conc * estc) + Math.Abs(cons * ests);
            goto _190;
            _180:
            result = conc * ress24 + cons * resc24;
            abserr = Math.Abs(conc * ests) + Math.Abs(cons * estc);
            _190:
            return result;
        }
    }
}